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Abstract 

We present a new first-principle theory for the calculation of the macroscopic second-order sus- 
ceptibility x i based on the Time-Dependent Density- Functional Theory approach. Our method 
allows to include straightforwardly the many-body effects. We apply the theory to the computation 
of the Second-Harmonic Generation spectroscopy, showing a very good agreement with experiment 
for cubic semiconductor GaAs. 
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Nonlinear optics is one of the most important and exciting field of fundamental and 
applied research, with applications in physics, chemistry and biology. Among the nonlin- 
ear phenomena existing in nature, the main role is played by second-order processes, like 
Second-Harmonic Generation (SHG) and Sum- Frequency Generation (SFG), which are ex- 
tremely versatile tools for studies of many kinds of surfaces and interfaces^. Their interest is 
rapidly growing, because of their exceptional sensitivity to space symmetry violations, and 
nowadays SHG and SFG are also used for characterizing systems like nanocrystal interfaces^ 
or as a probe for molecular chirality in polymer- and nanotubes^. The inverse process of 
SFG: Parametric Amplification together with Optical Rectification are used in microwave 
and terahertz technology^. 

In all these processes, the interaction of matter with light is described by the macroscopic 
second-order susceptibility x ■ 111 principle, x includes the many-body interactions 
among the electrons of the system: the variation of the screening fields on the microscopic 
scale, i.e. crystal local-field effects^ and the electron-hole interaction, i.e. excitonic effects^, 
as real and/or virtual excitations are created in the process. The theoretical description 
of the many-body effects in the second-order response is a big challenge and only a small 
number of ab initio works exist on the topic, mainly focused on the static or low energy 
limit. Self-consistent local-field effects were included in x^ within local density approxi- 
mation (LDA) through the "2n+l" theorem applied to the action functional as defined in 
Time-Dependent Density- Functional Theory (TDDFT)^ or in a band theory*^ in the case of 
semiconductor and insulator materials. Electron-hole interaction has been described, in the 
second-order response, through the solution of an effective two-particle Hamiltonian, derived 
in the Bethe-Salpeter equation approach (BSE)^ 1 ^. However, the validity of this approach 
has been demonstrated in linear-response calculations and the question arises whether it 
is possible to use this method also for higher-order calculations. This indicates the crucial 
need for different many-body approaches. 

In this letter we present a new first-principles theory for the calculation of the static and dy- 
namic macroscopic second-order susceptibility x i based on the Time-Dependent Density- 
Functional Theory (TDDFT) approach^ 1 ^. This formulation is valid for any kind of crystals: 
semiconductors and metals. The goal of our formalism is twofold: first, the exact relation 
(non-relativistic regime) between microscopic and macroscopic formulation of the second- 
order response, shown here for the first time, and second, a rigorous and straightforward 
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treatment of the many-body effects. In order to validate our theory we have applied our 
method to the SHG spectroscopy for the frequently studied material: cubic semiconductor 
GaAs. Indeed, GaAs has been object of large interest in SHG since its discovery, and a cer- 
tain number of experimental and theoretical studies have been performed, but no existing 
theoretical approaches have been able to give a conclusive comparison with experiments. We 
are able to do it with our formalism, showing here for the first time, an excellent agreement 
with experimental data. 

The objective of our theory is to find an expression for the susceptibility x i defined through 
the macroscopic second-order polarization 

E tot being the macroscopic component of the total electric-field. E tot includes the con- 
tribution from the external perturbing electric-field and from the electric-field due to the 
polarization of the system, induced by the external perturbation. 

Our starting point is the calculation of the second-order microscopic polarization via the 
second-order time- dependent perturbation theory. We obtain an expression for the micro- 
scopic polarization in term of the perturbing electric-field E p 

P (2) (q+ G,u) = J duj idtu 2 5(q. - qi - q 2 ) 

qi,q 2 ,Gi,G 2 

x5(u — Ui — w 2 )« (2) (q + G,qi + Gi,q 2 + G 2 ,Wi,w 2 ) 

xF( qi + G 1; Wl )£ p (q 2 + G 2 , oj 2 ), (2) 

where is the quadratic quasi-polarizability tensor. All quantities are functions of the 
frequency u, of the vectors q in the Brillouin zone and of the reciprocal lattice vectors G. 
In the first step of our calculation, all the quantities are written in term of the perturbing 
electric-field and are microscopic. At this point, there are two main difficulties in our 
formalism: first we need to average correctly in space to obtain macroscopic measurable 
quantities and second we have to express the polarization as a function of the total-electric 
field. 

To overcome these problems, we introduce the function 

a (1) (q, q,w) 



F(q,u) 



1 + 47T: 



47ra( 1 )' LL (q, q, u) 



(3) 



were oS 1 ^ is the linear quasi-polarizability and a^' LL is its longitudinal-longitudinal con- 
traction, and we use, in Eq. (J2]) the relation between the perturbing and the total electric 
field (macroscopic component) 

E^q,u) = F{q,u)E to \q,uj), (4) 

demonstrated by Del Sole and Fiorino^. 

We thus obtain the desired results: the macroscopic second-order polarization 



P$(q,u) = ^2 duidu 2 5(q - qi - q 2 ) 
qi,q2, 



x5(u -Ui- u 2 )F(q)a 2 (q, q ls q 2 , u l7 0J 2 ) 

xF( qi )F(q 2 )E tot (q u u 2 )E tot (q 2 ,u 2 ), (5) 

from which it is easy to derive, through a comparison with Eq. (pQ), our key quantity, the 
susceptibility ■ Pm (-^l- ©) contains the ah initio relation between the microscopic 
and macroscopic formulation of the second-order response, and it is shown here for the first 
time. 

Furthermore, our formalism is completely general for electric fields containing both longitu- 
dinal and transverse components. In the following, we will consider only the case of vanishing 
light wave vector (q — > 0), for which longitudinal and transverse responses are equal, since 
the direction of q is no longer defined for a uniform field and the responses depend only on 
the polarization of the fielcU^. Therefore, we have expressed the second-order susceptibility 
in terms of longitudinal quantities only. 

We have derived the optical susceptibility for any crystal symmetry and here we show the 
case of the cubic zinc-blend symmetry which has only one independent non- vanishing com- 
ponent (out of 18 tensor components) 

XxLi^i + wa, wi, u 2 ) = -^r lim 



xyzK .. x , ,•_,«, .oq r q ?/ q : 

xXp PP (q,q,^i,^2)eM (^i + u 2 )e L h f {uJi)e\f (u 2 ), (6) 

where q = (q x , q y , q z ). In this scalar equation, ej£ is the macroscopic longitudinal dielectric 
function, defined as D = e\jE tot (for longitudinal fields) and it has to be evaluated at the 
photon frequencies tui, u 2 and oj\ + uj 2 . 



The second-order susceptibility x X yz depends also on the xpppi which is the longitudinal 
second-order response function. All these quantities are evaluated in the limit q — > 0. The 
general expression for the Xpppi within TDDFT, is given through the generalized matrix 
Dyson equation 

Xo\u)i,U 2 ) [1 + fuxcMx (1) M] 

X [l + fuxcMX^M] 

+ Xo\"i + U2)9xc(ui + oo 2 )x {1) Mx {1) M (7) 

where f uxc is the sum of the bare-coulomb potential u and of the exchange-correlation 
kernel f xc = A new kernel, g xc = j^t, appears, defined as the second derivative of the 
exchange-correlation potential. Moreover, x (^0 i s the linear response function calculated 
via the Dyson equation 

[l-xS 1) («)/™c(«)]x (1) (w) = xS 1) (w)- (8) 

The functions xi^ i u ) an d X(P ( u ) are the linear and second-order response functions in the 
Independent-Particle Approximation (IPA). The response function Xppp can De calculated 
with different levels of approximation, depending on the kernels we use. Up to now, most of 
the ab initio calculations existing in literature were obtained within IPA, which we recover 
by setting f xc = and g xc = 0. In this case the factors 1 + ux^ and 1 — x<P u °f Eq. (JTj) 
compensate the ej^ functions of Eq. ([6]), leading to the usual expression Xxyz = Xo ■ Note 
that in Eq. (JTJ) and Eq. (jSJ) we have omitted the explicit dependence on the q and G-vectors 
of the response functions. 

To corroborate our theory, we show SHG spectra for cubic semiconductor GaAs, which is one 
of the most studied system in nonlinear optics from theory and experiments. We compare 
with the most accurate experimental data available now, obtained by Bergfeld and Daum^, 
who determined the \xxyz\ ( reported in Fig. (pQ)) over a wide spectral range. The abso- 

(2) 

lute value of the susceptibility \Xxyz\ has been determined by measurements of the reflected 
second-harmonic light (note that the presence of the surface leads to additional peaks that 
cannot be separated from the bulk contribution E=1.65 eV). In their work—, two theoretical 
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studies are merit ioned^ 1 ^, showing only a qualitative agreement with experiment, and not 
in the whole frequency range. Furthermore, none of them is able to reproduce the absolute 
values, even if many-body effects, in the framework of the BSE, are considered^. 
We have performed our electronic ground-state calculations with the plane-waves pseudopo- 
tential code ABINIT— and employed, for the nonlinear optics, a new nonlinear-response 
code implemented by us, on the basis of the linear-response Dp code2£. One of the crucial 
steps for the estimation of the second-order response function is the computation of the band 
structure of the system. For Gallium, Dal Corso et al.— demonstrate the importance of the 

(2) 

3d semi core states, pointing out that to obtain correct \Xxye\ values it is necessary to include 
at least the nonlinear core correction to the Gallium pseudopotential. In our calculations 
we have explicitly included the 3d semicore states of Gallium as valence states. In fact, the 

(2) 

presence or the absence of the 3d semicore states influences the value of \Xxyz\ for GaAs, 

(2) 

resulting in a slightly increasing of the \Xxyz\ when including the d states. 
In Fig. (Op), together with the experimental result^, is also shown the calculated \x£y \\ in the 
independent-particle approximation and with the inclusion of the crystal local-field effects. 
In both cases we have correctly^ applied the scissor shift^ of 0.8 eV to the Kohn-Sham band 
structure. Already within IPA we obtain the same shape as the experimental spectra. Both 
the theoretical and the experimental spectrum have three peaks and three valleys which are 
almost at the same energy position. This agreement can be achieved only when including 
explicitly the semicore states in the calculation. Taking into account the crystal local fields, 
for this energy range, the main effect is a decrease of the magnitude of the second-order 
susceptibility of the order of 10%, when compared to IPA. The shape of the spectrum is not 
significantly changed. However, even though the shape of the theoretical spectrum is good, 
the relative intensity of the peaks and in particular the magnitude of the susceptibility is 
not in agreement with the experimental values. The physics of the process is not sufficiently 
described neither in the independent-particle picture nor taking into account the microscopic 
inhomogeneity of the system. 

To go beyond these approximations, we have considered the excitonic effects through the f xc 
kernel, keeping g xc = 0. In Eq. (jSJ), f xc appears in the calculation of e^, an d Xppp- When 
using the ALDA kernel (not shown here) for f xc the result remains very close to those of IPA 
and IPA with crystal local fields. This is very similar to the failure of TDLDA absorption 
spectra in solids^, related to the lack of long-range contribution in the ALDA kernel. To 



6 



solve this issue, a model long-range kernel of the form a/q 2 has been proposed 2 ^, where a is 
a mean value for the dynamical dependence of f xc , in a given range of frequency. For GaAs 
the standard values for a are 0.05 in the static limit and 0.2 in the energy range of Fig. (J2J). 
The main effect of the a/q 2 kernel is to increase the magnitude of | Xxyz I recovering the order 
of magnitude of the absolute value of the experimental second-order susceptibility without 
changing the position of the energy peaks and valleys, as shown in Fig. fl2]). This behavior can 
be understood by solving analytically Eq. (J7]) without local fields, showing that the increase 
from the {x^l to the |x (2) | is proportional to [1 + a/4ir(e^(uj) - 1)] 2 [1 + a/4ir(e^(2uj) - 1)]. 
Only excitonic effects can correctly describe the magnitude of SHG measurements in GaAs. 
They have to be included carefully and consistently with a second-order process. In fact, 
in Fig. §2§ it is also reported the spectrum calculated by Leitsmann et al.— where excitons 

(2) 

are described in the second-order susceptibility using the BSE approach. The \Xxyz\ is much 
lower than the experiments. We believe that in the nonlinear response regime, for finite 
frequencies, the crystal local fields and the excitons, are not correctly described by this ef- 
fective Hamiltonian derived in the BSE approach. 

Instead, in the static limit, these effects seem to be less important and this BSE-based 

(2) 

method is still valid. Indeed, we obtain 205 pm/V for \x X yA which is in agreement with 
Chang et al. who obtained 236.4 pm/V and both results are in a reasonable agreement with 
experimental data: 180±10 pm/V— and 166 pm/V— at frequency 0.117 eV and 172 pm/V 
at frequency 0.118 eV— . 

In conclusion, in this letter, we have presented a new first-principles formalism for the cal- 
culation of the static and dynamic second-order susceptibility x^ 2 \ 111 t ne Time-Dependent 
Density-Functional Theory framework, valid for any type of crystal (semiconductor and 
metal). Our theory permits to write down, for the first time, the exact relation (non- 
relativistic regime) between microscopic and macroscopic formulation of the second-order 
response function, together with a rigorous and straightforward treatment of the many-body 
effects. We have applied this formalism to the most accurate experimental result made on 
cubic semiconductor GaAs, giving, for the first time, conclusive theoretical explanation of 
the experimental data. We show that only carefully including excitonic effects, it is possi- 
ble to recover the magnitude of the experimental x ■ Finally, we want to emphasize the 
importance of our formalism. This theory represents crucial progress for the development 
of nonlinear optics, which is a very active and exciting field in many disciplines of research 
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and we are convinced that this work can open new ways to explore it. 
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the Centre de Calcul Recherche et Technologie (CCRT). 
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FIG. 1: I Xxyz\ calculated within IPA (solid line) and including crystal local fields (dashed line). 
Comparison with the most accurate experimental data available nowi£ (circle). The experimental 
energy position of peaks and valleys are reported. 
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FIG. 2: The experimental Ixiyzl f° r GaAs (circle)^ is compared to our calculation (solid line), 
which includes the excitonic effects within TDDFT through the a kernel and with the calculation 
of Leitsmann et al.— (dot-dashed line) where the excitons are included within BSE framework. 



10 



